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FOREWORD 

This  report  documents  another  of  the  perturbing  terms  appearing  in  the  force 
equations  for  satellite  motion  in  the  TERRA  system  of  satellite  geodesy  computer 
programs.  A detailed  algorithm  is  presented  for  the  gravitational  action  which  the 
oceanic  tide  bulge  exerts  on  a satellite.  Two  earlier  reports  described  the  necessary 
equations  for  the  air  tides  caused  by  sun  and  moon.  Also,  the  manuscript  has  been 
completed  for  a further  report  which  will  present  an  algorithm  for  the  perturbing 
force  caused  by  the  tidal  deformation  of  the  solid  earth. 

The  following  text  was  reviewed  and  approved  by  Mr.  R.  J.  Anderle,  Head 
Astronautics  and  Geodesy  Division. 
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INTRODUCTION 


TERRA1  is  a new  system  of  computer  programs  for  satellite  geodesy.2  Like  its 
forerunners,  it  performs  geodetic  solutions  for  a variety  of  parameters  such  as  the 
constants  of  the  gravity  field  model,  the  observing  station  coordinates,  the  orbital 
parameters  of  the  satellites  involved  and  certain  of  the  geophysical  and  astronomical 
constants  which  determine  the  satellite  orbits.  To  be  able  to  take  advantage  of  the 
recent  progress  in  instrumentation  for  data  acquisition,  TERRA  is  a collection  of 


algorithms  some  of  which  are  rather  more  complex  than  the  corresponding  schemes 
from  its  forerunners.  In  particular,  the  force  equations  for  satellite  motion  now 
contain  terms  which  quite  realistically  indicate  the  effect  of  certain  physical 
phenomena  which  are  formerly  modelled  in  a crude  fashion  only  or  which  were 
discounted  altogether. 

Amongst  the  forces  previously  omitted  but  included  in  TERRA,  those  stand  out 
which  reflect  the  gravitational  action  of  the  various  types  of  earth  tide  on  the 
satellite  orbits.  They  are  the  air  tides  caused  by  the  moon  and  the  sun,  the  ocean 
tide,  and  the  tide  of  the  solid  earth  body.  Equations  for  the  two  atmospheric  tides 
were  elaborated  on  recently.3 ,4  Also,  the  solid  earth  tide  is  already  a part  of  the 
TERRA  coding.  It  is  the  subject  of  a technical  report.5 

For  the  solid  earth  tide,  we  obtained  a potential  from  the  literature  in  a form 
immediately  applicable  to  our  purpose.  Only  the  gradient  needed  to  be  found.  That 
involved  a large  computing  effort  which  was  in  the  end  accomplished  with  the  help 
of  an  analytical  computer  language.  For  the  two  air  tides,  surface  pressure  functions 
corresponding  to  the  atmospheric  tide  bulges  were  readily  available  in  the 
geophysical  literature.  From  surface  pressure,  we  managed  to  derive  the  disturbing 
acceleration  for  each  tide  via  Poisson  integration  followed  by  spatial  differentiation. 

The  present  report  on  the  ocean  tide  in  TERRA  is  the  third  in  a series  of 
four,  each  presenting  one  of  the  four  individual  tides.  Actually,  the  ocean  tide  was 
the  last  for  which  we  succeeded  in  our  effort  to  produce  an  adequate,  yet 
manageable,  disturbing  term.  For  this  there  are  reasons  inherent  in  the  problem. 
There  appears  to  be  a particular  aspect  to  the  ocean  tides  which  sets  them  apart 
from  the  remaining  earth  tides.  To  be  specific,  any  effort  to  mathematically 
formulate  the  shape  of  the  tidal  bulge  for  any  of  the  four  tides  may  be  expected 
to  involve  some  kind  of  eigenfunction  expansion.  The  same  is  true  for  any  de tailed 
description  of  items  derived  from  the  tidal  bulge,  such  as  Newtonian  potentials 
arising  in  connection  with  the  tidal  mass  redistribution  and  the  associated  gradients. 
For  the  solid  earth  tides  and  the  air  tides,  these  expansions  are  meaningful  on  a 
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domain  which  is  the  entire  earth  surface  and/or  the  space  above  it.  The  dominant 
boundary  condition  is  simple,  merely  requiring  continuity  and  periodicity  in 
longitude.  This  reflects  the  fact  that  there  exists  no  sharply  pronounced  barrier  to 
tidal  wave  propagation  in  the  atmosphere;  neither  is  there  such  a barrier  within  the 
solid  earth.  But,  barriers  of  this  type  do  occur  in  the  ocean.  They  are  the 
coastlines.  That  fact  gives  rise  to  a good  deal  of  mathematical  unpleasantness.  For 
example,  if  one  desires  to  represent  the  tidal  ocean  surface  by  the  customary 
surface  harmonic  expansion,  one  will  have  to  decide  whether  the  mathematical 
expression  for  that  surface  is  to  be  valid  for  any  possible  pair  of  values  for  latitude 
and  longitude  or  if  one  wishes  to  produce  an  expression  which  is  valid  only  for  a 
domain  which  coincides  with  the  ocean.  The  former  expression  will  inherently  be 
capable  of  yielding  zero  tidal  deflection  for  any  point  located  inland,  within  a 
specified  threshold.  But  it  obviously  will  have  to  contain  a formidable  number  of 
parameters.  The  latter  expansion  will  be  of  the  type  which  the  physicist  calls 
“non-analytic.”  It  will  be  simpler  than  the  first,  but  it  will  be  valid  only  on  the 
ocean  domain.  Thus,  we  shall  have  to  augment  it  by  the  familiar  postulate  that  it 
be  disregarded  on  that  domain  which  consists  of  land  and  that  its  value  there  be 
zero  by  definition. 

We  did  indeed  encounter  the  latter  situation  when  we  tried  to  adapt 
Hendershott’s  expansion  for  the  ocean  tide  surface6  to  our  problem.  On  a superficial 
inspection,  Hendershott’s  equation  for  the  tidal  surface  appeared  well-suited  for 
conversion  to  a tidal  mass  layer,  as  done  successfully  for  the  air  tides  in 
References  3 and  4.  But,  when  we  subsequently  tried  to  actually  formulate  and 
execute  the  Poisson  integral  to  calculate  the  tide  potential  exterior  to  the  earth  as 
we  had  done  for  the  air  tides,  we  immediately  faced  the  task  of  spelling  out  the 
integration  boundary.  This  task  looked  difficult  because  Hendershott's  -equation 
appeared  to  be  of  what  we  called  the  “non-analytical”  type  above.  Consequently,  we 
had  the  option  of  trying  to  develop  integration  boundaries  which  would  reflect  the 
continental  coastlines.  Or  else,  we  might  attempt  to  achieve  very  simple  integration 
boundaries  by  converting  the  formula  for  the  tidal  surface  into  one  which  would  be 
valid  on  the  entire  globe,  yielding  spurious  values  of  negligible  magnitude  on  land. 
Because  of  the  great  complexity  inherent  in  either  approach,  we  abandoned  the 
Hendershott  model.  This  decision  was  made  easy  by  the  circumstance  that  the 
number  of  expansion  coefficients  listed  in  Reference  6 was  clearly  insufficient  for 
our  purpose.  Further,  our  reference  indicated  that  an  additional  number  of  these 
coefficients  had  actually  been  computed.  But  we  were  unable  to  obtain  those  from 
the  author. 
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In  fact,  we  abandoned  not  only  the  Hendershott  model,  but  we  resolved  to 
entirely  stay  clear  of  concepts  like  continuous  surface  densities  and  Poisson  integrals. 
We  were,  of  course,  aware  of  the  important  work  done  in  the  field  of  earth  tide 
effects  on  satellite  motion  which  is  being  published  under  the  names  of  K. 
Lambeck,  A.  Cazenave,  and  G.  Balmino.7,8,9  We  hesitated  to  get  involved  with 
these  papers  because  they  represent  continuing  research.  They  certainly  identify 
various  significant  components  of  the  tide  potential.  We  desired,  however,  to  specify, 
for  use  with  TERRA,  a tide  potential  which  would  implicitly  contain  not  only  the 
known  tidal  effects  but,  additionally,  as  many  of  those  perturbations  which  might 
be  characteristic  of  any  particular  type  of  satellite  orbit  to  which  TERRA  might  be 
applied. 

Fortunately,  an  entirely  different  approach  offered  itself.  A new  theory  for  the 
ocean  tide  was  recently  developed  at  the  Naval  Surface  Weapons  Center,  Dahlgren 
Laboratory  (NSWC/DL),  by  Dr.  E.  Schwiderski.  This  is  an  improved  Zahel  theory.* 
A computer  model  which  implements  this  theory  has  just  begun  to  yield  results 
which  are  in  excellent  agreement  with  the  observed  data.  At  the  present  time,  this 
model  reflects  only  the  M2  component  of  the  tide.  The  latter  is  thought  to  account 
for  a substantial  portion  (about  70  percent)  of  the  oceanic  mass  dislocation  due  to 
the  tides.  Also,  it  appears  that  it  is  not  a technical  problem  but  purely  a question 
of  funding  to  generalize  this  model  to  include  any  desired  member  of  the  tide 
spectrum  for  which  sufficient  observations  are  available.  As  the  present  version 
appears  by  itself  to  be  the  most  accurate  and  detailed  one  amongst  the  quantitative 
descriptions  known  for  the  ocean  tide,  we  confidently  based  our  work  on  it. 

The  Schwiderski  model  presents  the  tidal  ocean  surface  as  a listing  of  values, 
for  the  tidal  amplitudes  and  phase  angles,  at  the  center  points  of  surface  area 
elements  which  form  a grid  covering  most  of  the  oceans.  This  permitted  us  to 
specify  two  algorithms  which,  if  used  in  sequence,  should  enable  TERRA  to  account 
for  the  orbital  perturbations  due  to  the  tide  bulge. 

In  detail,  we  regard  the  tidal  height  on  each  ocean  surface  element,  as 
computed  by  the  Schwiderski  model,  as  a measure  for  that  portion  of  the  mass  of 
the  tidal  bulge  which  is  located  within  the  particular  area  element.  We  postulate 
that  each  tidal  height  value  thus  resembles  a gravitating  point  mass  which  perturbs 
the  satellite  by  its  presence.  The  gravitational  potential  for  each  point  mass  is  then 
expanded  in  terms  of  Legendre  functions.  Upon  further  summation,  over  all  area 
elements,  the  resulting  perturbing  potential  is  decomposed  into  products  of  terms 
which  contain  the  source  point  coordinates  and  other  terms  which  are  functions  of 


*For  details  on  the  Zahel  theory,  see  References  10  and  II. 
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the  field  point  (satellite  position)  and  time  only.  The  latter  terms  act  the  part  of 
eigenfunctions  while  the  former  are  the  expansion  coefficients  of  the  perturbing 
potential. 

The  quantity  actually  needed  is  the  perturbing  acceleration  which  is  the 
gradient  of  the  perturbing  potential.  This  can  be  easily  assembled  from  the  just 
mentioned  expansion  coefficients  and  eigenfunctions  because,  as  any  reasonably 
complete  text  on  potential  theory  or  related  subjects  will  show,  it  is  possible  to 
express  the  spatial  derivatives  of  the  solid  harmonics  occurring  in  our  potential  as 
linear  combinations  of  certain  of  the  solid  harmonics  themselves. 

The  resulting  scheme  for  computing  the  perturbing  acceleration  quite  closely 
resembles  that  used  to  evaluate  the  spherical  harmonics  occurring  in  that  part  of  the 
TERRA  equations  of  motion  which  deals  with  the  main  part  of  the  earth’s 
gravitational  field.  Every  effort  was  made  while  constructing  the  two  algorithms  to 
enable  the  programmer  to  utilize  relevant  procedures  for  which  codings  already  exist 
such  as  the  recurrence  relationships  for  the  spherical  harmonics.  The  latter  are  used 
in  a number  of  computer  programs  for  satellite  geodesy1 2 because  they  greatly 
facilitate  the  computation  of  the  eigenfunctions  occurring  in  the  expressions  for  the 
geodetic  potentials,  and  thus,  in  the  perturbing  accelerations  which  act  on  the 
satellites  as  a consequence  of  the  presence  of  these  potentials. 

What  might  conveniently  have  been  one  single  collection  of  equations  for 
coding  was  partitioned  into  two  separate  algorithms,  for  various  reasons.  Work  done 
in  the  past1 2 had  demonstrated  the  power  and  convenience  of  our  recursive  schemes 
for  the  eigenfunctions  of  the  potential.  These  codings  require  that  the  constituents 
of  the  potential,  namely  the  eigenfunctions  and  the  expansion  coefficients,  be  known 
before  the  field  gradient  may  be  evaluated.  Finding  the  expansion  coefficients  and 
the  eigenfunction  values  is  the  tedious  part  of  the  task.  But,  once  they  are  known, 
the  the  corresponding  gradient  vector  follows  easily.  This  suggested  that  we  relegate 
the  two  just  mentioned  operations  to  a preprocessor  algorithm  which  feeds  data  into 
a second  algorithm  which  is  part  of  the  main  body  of  TERRA  and  which 
constitutes  the  perturbing  acceleration  associated  with  the  tide.  Also,  a second 
motive  for  dividing  the  procedure  was  that  we  anticipated  that  we  might  wish  to 
study  certain  aspects  of  the  tide  potential  which  arc  rather  unrelated  to  satellite 
motion  in  TERRA.  One  obvious  application  would  be  to  find  out  how  well  the 
original  tide  bulge  can  be  resynthesized  from  the  potential  generated  by  the  above 
mentioned  point  masses.  Investigations  of  this  type  would  be  greatly  aided  if  it  were 
possible  to  isolate  that  coding  which  produces  the  ocean  tide  perturbing  potential 
from  the  body  of  TERRA. 
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As  will  be  seen  below,  this  separation  was,  by  necessity,  not  a perfect  one. 
But,  it  provides  for  a preprocessor  which  contains  most  of  the  individual  features  of 
the  tide  model  in  use,  such  as  the  parameters  of  the  tidal  point  masses.  At  the 
same  time  it  permits  the  actual  perturbing  term  to  appear  in  the  equations  of 
motion  in  a form  which  is  quite  independent  of  the  individual  properties  of  the 
tide  model. 

More  precisely,  the  first  of  our  algorithms  stands  alone,  by  itself,  outside  the 
body  of  computer  programs  which  form  TERRA.  Its  task  is  to  find,  from  the  tidal 
distortion  of  the  ocean  surface,  the  numerical  values  of  the  tidal  point  masses  as 
well  as  their  positions.  From  here,  it  calculates  the  expansion  coefficients  for  the 
perturbing  potential.  Also,  it  is  capable  of  calculating  values  for  the  eigenfunctions 
of  the  potential,  for  any  given  field  point. 

Along  with  items  like  the  tidal  frequency  and  related  astronomical  data 
(specifying  the  position  of  the  celestial  body  or  bodies  causing  the  tide),  the  just 
mentioned  expansion  coefficients  are  the  input  for  the  second  routine.  The  latter  is 
an  integral  part  of  the  TERRA  equations  of  motion.  Like  the  first  algorithm,  it  is 
equipped  to  evaluate  the  eigenfunctions  (which  are  now  to  be  regarded 
“eigenfunctions  of  the  perturbing  potential”).  The  results  of  that  procedure  are 
subsequently  combined  with  the  expansion  coefficients  to  form  the  spatial 
components  of  the  perturbing  acceleration  for  the  satellite  motion.  Note  that  the 
separation  of  the  two  algorithms  is  imperfect  because  there  appear  certain  multipliers 
in  both  of  them  which  reflect  the  time-dependent  nature  of  the  tide.  In  the  form 
in  which  they  appear  below,  these  multipliers  are  valid  for  the  M,  tide  only. 
Caution  will  have  to  be  exercised  to  adapt  them  to  any  different  tide  which  may 
be  introdueed  later. 

The  first  algorithm  will  be  executed  infrequently.  As  the  Schwiderski  model 
specifies  its  results  in  terms  of  discrete  values  at  grid  points  which  are  1 degree 
apart  in  longitude  and  latitude,  our  first  algorithm  utilizes  about  45  000  tidal 
amplitudes  and  phases  to  compute  45  000  complex  numbers  (number  pairs) 
representing  point  masses  which  are,  in  turn,  converted  to  45  000  pairs  of  complex 
spherical  harmonic  coefficients.  One  or  two  dozen  of  these  pairs  of  complex 
spherical  harmonic  coefficients  will  be  hand-selected  for  input  to  CELEST/TERRA 
which  will  be  modified  to  operate  on  these  coefficients  with  the  second  algorithm. 
Experiments  will  then  be  conducted  to  determine  the  minimum  number  of 
coefficients  needed  to  reflect  the  action  of  the  ocean  tides  on  the  satellite  orbits. 
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TIDAL  POINT  MASSES 


To  find  the  tidal  point  masses  and  their  positions,  divide  that  portion  of  the 
globe  which  is  ocean  into  area  elements,  according  to  the  output  format  of  the 
Sehwiderski  Ocean  Tide  Program.  For  details,  see  Figures  1 and  2.  Note  that  each 
surface  area  element  has  sides  the  angular  length  of  which  is  1 degree.  Exclude 
from  consideration  all  those  area  elements  which  the  Ocean  Tide  Program  discounts, 
especially  those  located  over  land  and  all  those  which  are  located  below  a certain 
southern  latitude.  These  area  elements  which  contain  the  North  pole  are,  of  course, 
triangles.  For  the  purpose  of  mass  point  location,  assume  that  the  area  elements  be 
located  on  the  reference  ellipsoid.  As  far  as  area  size  is  concerned,  assume  that  the 
area  elements  be  located  on  a sphere  of  radius  R.  Most  frequently,  R will  be 
equated  to  the  semimajor  axis  of  the  reference  ellipsoid.  Specify  R in  terms  of 
kilometers. 

Now,  place  a point  mass,  nVj,  into  the  center  of  all  “valid”  area  elements,  as 
indicated  in  Figure  2.  Number  all  surface  area  elements,  point  masses  and  associated 
quantities  either  by  a double  subscript,  ij,  or  index,  v,  as  expedient.  Generally, 
regard  ij  and  v as  interchangeable  indices.  Figures  1,  2,  and  3 now  suggest  how  to 
form  the  expressions  for  the  surface  areas  of  the  area  elements, 


AS  = 


7T  \*  , / 7T 

R"  sin  ( j 

180/  V 180 


j = 2,3,4,  • ••  ,jmax  < 180, 


= 1 (-S*2 

2 \ 1 80  / 


which  is  the  surface  area  of  the  “polar”  surface  element. 
Specify  the  positions  of  the  point  masses  as  follows. 


i?  = j 

> 180 


Xjj  = (i  2) 


AREA 


AO 

2 

1. 


////////; 


ELEMENT  ij 


M 


i,J 


AX 

2 


'"/////////'///S'" 


\y/// /////// ///////', 

AX 


x = (i-l) 


X = Xj  = i 


Figure  2.  Position  of  Point  Mass  at  the  Geometrical  Center  of  the 
Surface  Area  Element,  ij 
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Figure  3.  Position  of  the  Point  Mass  in  the  Earth-Fixed  Cartesian 
Coordinate  Frame  y{1),  y(2>,  and  yl3>  [where  y(11  is  in  the 
Direction  of  the  Point  where  Greenwich  Meridian  and  Equator  Intersect  and 
y121  is  in  the  Equatorial  Plane] 
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the  colatitudes  associated  with  the  area  element,  ij. 


a 

c 


L 


o = o 

ij 


V 


colatitude  associated  with  the  point  mass,  ij. 

7T  7T  7T  / 1 \ 

0 . = - - 0.  = - - j - - 

■J  2 'J  2 180  V 2/ 


(geocentric)  latitude  of  m ^ . 


longitude  (east)  of  m^. 


(4) 


(5) 


(6) 


(7) 


the  geocentric  distance  of  m(j . 


According  to  Figure  3,  the  earth-fixed  Cartesian  coordinates  of  are 


xij  = Pjj  cos  0 j . cos  X(j 

(8) 

Yij  = Pij  cos  0jj  sin  Xjj 

(9) 

ZU  = sin  °ij 

(10) 

Above,  e2  is  the  square  of  eccentricity  of  the  reference  ellipsoid.  In  case  it  is 
desired  to  start  from  the  flattening,  f,  of  the  ellipsoid,  find  e2  from 


e2  = (2  - f)f 


(11) 
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'..js  , .it 


r,  • 


A * 


Further,  p is  the  density  of  sea  water  and  G is  Newton’s  constant  of  mass 
attraction.  Assume 


p = 1 


(12) 


for  a sea  water  density  of  1 g/cm3 * * * . If  one  wishes  to  employ  a different  density 
value,  scale  p accordingly.  Also, 


G = (6.67  E-  20)  km3  /kg  sec2 


(13) 


Now,  let  be  the  tidal  amplitude,  in  meters,  on  area  element,  ij.  6^  is  the 
associated  tidal  phase  angle.  Both  constitute  printout  from  the  Ocean  Tide  Program, 
t*  is  Universal  time,  in  seconds  (reference  zero  is  midnight  UT).  a is  the  circular 
frequency  associated  with  the  M2  component  of  the  ocean  tide  (two  cycles  per 
lunar  mean  day). 

X is  the  mean  mean  longitude  of  the  moon,  in  degrees,  at  the  beginning  of  the 
particular  day  UT.  A routine  is  indicated  for  x>  below.  And 

1 80 

a - (1.40519E-  04)  sec-1  (14) 

7 T 


Then,  define  the  “gravitational  charge”  (capable  of  assuming  positive  as  well  as 
negative  values)  of  the  point  mass: 

1 

m„  = — [«,  cos  (at*  + x)  + 0„  sin  (at*  + x)J  (15) 

Cj 

a,  = 109  p G AS^  cos  8v  (16) 

= 109  p G ASy  sin  dv  (17) 


expected  to  result  in  terms  of  km3/sec1 2.  Thus,  the 


1 1 


The  latter  two  quantities  are 
unit  for  will  be  kilograms. 


Let  N be  the  number  of  point  masses.  We  expect  N to  roughly  equal  45  000. 


To  find  x.  execute  the  following  routine.  Let  J be  the  number  of  the  calendar 
year  (1975  or  1976  or  1977;  generalization  to  other  years  is  obvious),  M be  the 

number  of  the  individual  month  in  the  calendar  year,  and  D be  the  number  of  the 

particular  day  in  the  month.  Also,  let 

1 n = m 

6(n,m)  = if  (18) 

0 n ¥=  m 

Then 


r 


max  n 

0 = V !'  (F  U + H V 

^ ,,  ,,  n m nm  n in  n in 

n = 0 in  =0 


(26) 


Fnm  and  Hnm  are  the  expansion  coefficients.  Unm  and  Vnm  are  the 
eigenfunctions.*  Note  that  the  coefficients  of  the  expansion  can  be  expected  to 
depend  on  time,  because  the  point  masses  from  which  they  result  oscillate  with  the 
tide.  The  eigenfunctions  are  solid  harmonics  which  in  turn  are  functions  of  the  field 
point  coordinates  (satellite  position  vector). 

Obtain  the  necessary  input  data  from  the  previous  chapter.  Note  that  p is  the 
earth’s  gravitational  constant.  Select  for  it  a value  compatible  with  the  rest  of  the 
TERRA  program.  Express  p in  terms  of  kilometers  and  seconds.  Find  now 


F = a,.  cos(at*  +x)  + 0F  sin(ot*  + x) 

II  III  * n m 1 n ni 


H. 


a.. 


aH  cos  (at*  + x) +I^h  sin(ot*+x) 

n n m n m 


■ (2 " 5") 


2 - 5"  I <n  ~ m>1  — ” 

(n  + m)!  pR2n  „ = i 


2 p2n+1  a r 

^ » v n m 


PF  = (2-  5"^  01  ■ m>:  ~ V p- n + 1 p fv 

nm  V "'/  (n  + ni )!  pR2n  1 v nr 


Sp2n  + ,a  h- 
(n  + m)!  pR2n  „=i  v v nr 


Ph  , 


(n-m)!  1 


- 2 p2n+l  3 hv 

n _ , ^ v n 


r 

n m 


hnm 


(n  + m)!  pR  „=i 

= -^Tp:(sin^)cosmX, 

^ V 
Rn 

= _ . - Pm(sin0  ) sin  mX 

„n  + 1 n v v 


5m  = 


1 n = m 

if 

0 n ¥=  m 


For  details  on  the  eigenfunctions,  see  the  next  chapter. 
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(27) 

(28) 

(29) 

(30) 

(31) 

(32) 

(33) 

(34) 

(35) 


Find  the  f"  and  h"  as  follows.  Noting  that 

n m n in 


h'.o  = 0 


for  all  values  of  n,  calculate  now,  separately  for  each  v,  the  required  f*m  and  h^t 
from  the  following  recurrence  relations.  To  advance  in  n,  evaluate 

fn%.  ,m  = — ^71  [(2n  + ‘ } Sil1  " ("  + >•'"]  <3> 

hn  + 1 ,m  = ^n+DsineX.™  -(n  + m)PX-iJ  (3< 

To  advance  m,  use 

fn  + I ,n  + 1 = (2n+  c°sXvr  n -COS0,  Si,,\hn.n]  (4< 

hn  + i.n  + i = (2n  + Dp,,  [cose,,  cos  X„h^  n + cos  0,  si"  MnJ  (4 

Start  these  recurrences  from 


hM  = 0 


■ 


- 


EIGENFUNCTIONS  FOR  THE  TIDE  POTENTIAL 


The  eigenfunctions  for  the  tidal  potential  are 


unm  = ^Ti  P^'  (sin  0)  cos  mX 
uRn 

Vnm  = ^T-1  P"1  (sin  0)  sin  mX 


(45) 

(46) 


Let  y(l),  y(2),  and  y(3)  be  the  earth-fixed  Cartesian  coordinates  (see  Figure  3)  of 
the  general  point  (“field  point”).  Note  that 


r 


= + 


(yO)2  + (y(2))2  +(y(3))2 


1/2 


(47) 


Also  observe 


R 


P = 


(48) 


, ym 

sin  i//  = 

r 


(49) 


and 


V 


n ,0 


= o 


(50) 


for  all  values  of  n.  Now  calculate  (numerically)  the  required  U and  V from 

* n in  n ni 

the  following  recurrence  relations. 

To  advance  in  n,  evaluate 


U 


n + 1 .m 


[(2n  + 1 ) sin  i^U 

n - m + 1 n,m 


- (n  + m)pUn_  , m 


(51) 


I 

e 


15 


hL 


r*i- 


> . V 


[(2n  + 1 ) sin  i^Vn.m  - (n  + m)pV 

n - m + 1 


(52) 


V 


n +■  I ,m 


n - 1 ,m 


1 


To  advance  m,  use 


"yd) 

v<2) 

^n+  1 ,n  + 1 

= (2n  + l)p 

r 

u - - — 

n ,n  r 

v, 

"y(D 

y(2) 

v 

v n + 1 ,n  + 1 

= (2n  + 1 )p 

r 

V + - — 

n ,n 

u 

1 

(53) 

(54) 


Start  from 


U, 


0,0 


U 


l ,o 


Ry(3) 

= M 


V = V 

0,0  vl,0 


(55) 

(56) 

(57) 


A TIDE  POTENTIAL  COMPUTER  ROUTINE 


This  is  the  first  of  the  two  algorithms  mentioned  in  the  introduction.  It  is 
intended  to  be  a multipurpose  computer  program,  capable  of  calculating  the 
Newtonian  potential  associated  with  the  ocean  tide.  We  propose  that  it  be  named 
Ocean  Tide  Potential  Routine.  This  routine  has  three  segments,  each  indicated  by 
one  of  the  preceding  three  chapters. 

The  first  program  segment  is  supposed  to  be  coded  from  Equations  1 through 
25.  As  already  explained,  it  converts  the  tidal  heights  into  point  masses,  complete 
with  expressions  for  the  “gravitational  charge”  and  associated  position  vector. 
Observe  that  the  m^  oscillate  with  the  tide  (Equation  15).  It  must  never  be 
forgotten  that  several  of  the  equations  apply  strictly  to  the  M,  (semidiurnal  lunar) 
tide.  Affected  are  those  equations  which  contain  a and  x-  Also,  of  course,  the  tidal 
amplitudes,  fjj,  and  the  5jj  are  valid  for  the  M2  tide  only.  Should  it  be  desired  to 
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calculate  the  disturbing  potential  associated  with  a different  ocean  tide,  a suitable- 
table  will  have  to  be  obtained  from  a theory  for  that  tide  for  the  and  6(j.  Also, 

the  new  tide  frequency  and  a new  algorithm  for  the  position  of  the  tide  generating 
celestial  body  will  have  to  be  substituted  for  our  present  a and  x- 

The  second  program  segment  is  to  be  assembled  from  Equations  26  through  44. 
It  calculates  the  expansion  parameters  for  the  potential.  Note  that  this  involves 
certain  vast  summations  over  the  individual  point  masses.  Independent  of  the  nature 
of  the  application  contemplated  for  the  Tide  Potential  Routine,  these  summations 
must  be  extended  over  all  surface  area  elements.  Neither  physical  nor  mathematical 
criteria  exist  which  would  permit  us  to  make  any  exceptions.  The  number  N of 
mass  points  is  equal  to  the  number  of  area  elements  (each  of  the  latter  being 
1 degree  square  in  size)  necessary  to  cover  the  globe  except  for  the  land  areas  and 
also  excluding  a substantial  cap  centered  upon  the  South  pole  (as  specified  by  the 
Ocean  Tide  Program).  N is  about  equal  to  45  000.  However,  only  four  summations 
are  involved  and  these  need  to  be  done  only  once,  because  neither  the  tide 
amplitudes  nor  the  phases  depend  on  time. 

The  third  algorithm  segment  consists  of  Equations  45  through  57.  This  is  for 
the  calculation  of  the  eigenfunctions.  To  avoid  a misunderstanding,  we  emphasize 
that  it  is  strictly  a numerical  scheme.  No  analytical  expressions  are  being  generated 
for  the  eigenfunctions.  Rather,  this  portion  of  the  Tide  Potential  Routine  starts 
from  values  for  the  field  point  coordinates  and  proceeds  by  numerically  evaluating 
the  individual  eigenfunctions.  Because  this  is  done  recursively,  all  eigenfunctions  must 
actually  be  evaluated  for  any  particular  position  vector,  up  to  a limit  n on  the 
order  and  the  rank  of  the  functions,  irrespective  of  the  nature  of  the  application. 
Note  that  n is  the  upper  summation  limit  in  Equation  26.  It  depends,  of 
course,  on  the  particular  application. 


PERTURBING  ACCELERATION 


The  perturbing  acceleration  to  be  inserted  into  the  TERRA  equations  of  motion 
is  the  positive  gradient  of  the  ocean  tide  potential.  Execute  now  the  following 
procedure,  separately  for  each  time  line  in  the  orbit  integration. 

For  each  time  line,  x(l),  x,2).  and  x13’  are  the  Cartesian  components  of  the 
satellite  position  vector  in  inertial  space.  Assume  that  the  x("  are  given  in  terms  of 
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kilometers.  Part  of  the  TERRA/CELEST  equations  of  motion  is  an^  algorithm  which 
performs  the  necessary  transformations*  between  the  inertial  frame,  x"\  and  the 
earth-fixed  frame,  y<u.  Now  apply  this  algorithm  to  rotate  the  inertial  position 
vector,  x(i\  into  the  corresponding  earth-fixed  position  vector,  ytl).  Eind 

r = +{(xtn)2  +(xl2))2  +(x(3>)2}  1/2  (58) 


or 


r = + |(yl  1 *)2  + (y(2))2  + (y<3))2 } 1/2  (47) 

as  convenient.  Calculate  the  values,  associated  with  the  satellite  position  vector,  y"\ 
of  the  eigenfunctions  Unm  and  Vnm. 

Evaluate  now  the  Cartesian  components  of  the  perturbing  acceleration,  in  the 
earth-fixed  frame: 


dip 

nm  ax 
2 

i ( 

t 

F 

dUnm 

+ H , 

9Vnni\ 

(59) 

9y(1> 

n = 0 

m = 0 

y n m 

9y(1) 

nm 

9y(1)  / 

dip 

9 y ( 2 ) 

"max 

2 

n = 0 

n 

I 

m = 0 

(F„,» 

aunm 

9y(l  ) 

+ Hnm 

9Vnm\ 
9y(2)  / 

(60) 

dip 

n ni  a x 
2 

n 

2 I 

fF 

9Unm 

+ H 

dV„m\ 

(61) 

9y(3  ) 

n = 0 

m = 0 

i 1 nm 

9y(  1 ) 

n m 

9yO)  J 

where 


+ 1 ,m  + 1 


‘References  for  this  transformation  are  unavailable  at  present.  But  the  transformation  is  already  part  of  the  IfRRA 
coding  and  will  appear  in  the  formal  TERRA  documentation. 

**That  is  the  frame  in  which  the  orbit  integration  is  performed. 


dUnm 

9y(  I ) 


1 M f _ l 

l -)  *\n  n ^n  + 1 ,m  - 1 2 
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1 / 1 1 

— I-  -A  V - -V 

I ■)  mn  n + l,m-l  *>  n + l,ni+l 


- - (n  - m + 1)U  .. 

' n + I ,m 


1 \ 1 
— -A  V - - v 

7 in  n n + 1 ,m  - 1 9 n + 1 ,m  + I 


1/] 

l 7 n + 1 ,m  - 1 r -)  + 1 ,in  + 1 


(n  - m + 1 )Vn  + , „ 


and  where 


Amn  = (n  - m + l)(n-  m + 2) 

u = - (n  - 1)!  U 
(n  + 1 )! 

(n-  !)! 

n-1  (n  + 1 )!  "-1 


To  find  the  Cartesian  components  of  the  perturbing  acceleration  in  the  inertial 
frame,  rotate  the  vector  3$/3y(l)  (i  = 1,  2,  3)  back  into  the  x(l)  frame. 


OCEAN  TIDES  IN  THE  TERRA  EQUATIONS  OF  MOTION 


The  Equations  58  through  70  listed  in  the  preceding  chapter  constitute  the 
second  of  the  two  algorithms  mentioned  in  the  introduction.  It  might  suitably  be 
named  Ocean  Tide  Disturbing  Acceleration. 


Formula  26  for  the  tide  potential  involves  a summation  over  the  index  n,  to 
be  extended  to  the  summation  limit  n . In  case  the  Ocean  Tide  Potential 
Routine  is  employed  to  conduct  the  previously  mentioned  experiment  of 
resynthesizing  the  tidal  bulge,  a large  number  of  terms  will  have  to  be  included  in 
the  potential.  To  be  specific,  the  number  of  terms  will  have  to  be  about  equal  to 
the  number  of  point  masses.  The  latter  number  is  N.  n x can  be  expected  to  be 
a number  near  to  the  square  root  of  N.  However,  when  applying  the  Ocean  Tide 
Potential  Routine  to  the  computation  of  the  perturbing  acceleration  due  to  the 
ocean  tide,  in  THRRA,  a comparatively  small  number  of  terms  should  be  sufficient. 
Only  those  terms  should  then  be  included  in  the  potential  which  are  significant  for 
the  purpose  of  the  orbit  computation  (short  arcs  and  long  arcs  of  the  satellite 
trajectory  each  suffer  perturbations  which  may  arise  from  different  members  of  the 
perturbing  potential  expansion).  We  shall  not  try  to  identify  the  useful  terms  now, 
a priori.  Instead,  we  have  chosen  to  proceed  empirically.  As  far  as  their  application 
to  THRRA  is  concerned,  we  are  hereby  submitting  both  algorithms  for  coding,  as 
they  stand,  preserving  full  generality  as  far  as  the  retention  or  deletion  of  any  terms 
is  concerned.  We  are  proposing  that  this  latter  question  be  resolved,  later  on,  by 
exploratory  THRRA  or  CHLEST  runs. 

The  reader  may  have  noticed  that,  directly  following  Equation  47  in  the 
preceding  chapter,  we  specified  that  the  eigenfunctions  be  evaluated  for  the 

earth-fixed  satellite  position.  It  was  implied  that  this  should  be  done  by  using  the 
recursive  scheme  for  the  eigenfunctions  outlined  in  the  Ocean  Tide  Potential 

Routine.  We  now  propose  that  the  programmer  decide  whether  he  wishes  to  call 
upon  the  Ocean  Tide  Potential  Routine  every  time  its  use  is  needed  during  the 
execution  of  the  Ocean  Tide  Disturbing  Acceleration  routine.  Remember  that  the 
Ocean  Tide  Potential  Routine  is  external  to  TERRA  while  the  Ocean  Tide 
Disturbing  Acceleration  is  part  of  the  THRRA  coding.  We  imagine  that  the 

programmer  may  prefer  to  duplicate  the  scheme  for  calculating  the  eigenfunctions 
and  make  it  a part  of  the  coding  of  the  Ocean  Tide  Disturbing  Acceleration.  That 
would  make  the  ocean  tide  term  an  entirely  self-contained  part  of  THRRA. 
analytically.  Of  course,  the  expansion  coefficients  would  still  have  to  be  obtained 
from  the  Ocean  Tide  Potential  Routine  preprocessor. 

Finally,  one  more  item  remains  to  be  discussed  which  is  related  to 

computational  economy.  This  concerns  the  time  dependent  factors 


(at*  + x) 


• Vv . 'r  .i$*h>  ..-»**  . ».  f 


occurring  in  the  expansion  coefficients  for  the  tide  potential.  Because  these  terms 
are  present  in  the  tide  potential,  they  will  also  appear  in  the  perturbing  acceleration 
associated  with  the  ocean  tide.  Consequently,  it  will  be  necessary  to  cope  with  them 
during  each  step  of  the  (numerical)  orbit  integration.  In  turn,  the  frequent  need  to 
evaluate  these  trigonometric  expressions  can  be  expected  to  very  adversely  affect  the 
time  required  to  perform  the  orbit  integration.  We  thus  propose  that  it  be  avoided 
to  calculate  these  factors  separately  for  each  time  line.  Instead,  the  following  scheme 
may  be  followed. 

Let  t*  and  t* , be  the  values  of  Universal  time  for  which  subsequent 

i \+ 1 

integration  steps  (time  lines)  are  to  be  performed.  Update  the  trigonometric  time 
factors  as  follows. 


cos  (at*  , + X)  = COS  | (at*  + x)  + «At*  J = cos  (at*  + x)  cos  aAtH 
-sin  (at*  + x)  sin  aAt* 


sin  (at*  , + x) 


sin  (at*  + x)  cos  aAt*  + cos  (at*  + x)  sin  aAt* 


At*  = t*  , - t* 


Using  the  algorithm  for  x (Equations  18  through  25),  update  x whenever  t*  enters 
the  following  day. 
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APPENDIX  A 

NOTES  ON  THE  TIDE  POTENTIAL 


Two  additional  comments  may  be  useful  in  connection  with  the  tide  potential. 
The  first  concerns  some  details  of  the  derivation.  The  second  deals  with  the 
frequently  mentioned  recurrence  relations  for  the  expansion  coefficients  and 
eigenfunctions. 

When  deriving  the  tide  potential,  we  started  from  the  potential  function  for  the 
individual  point  mass. 


(74) 


where  r is  the  position  vector  of  the  general  point;  pj;  indicates  the  position  of  the 
point  mass.  Let 


(75) 


Let  r,  0 , and  X be  the  geocentric  distance,  geocentric  latitude,  and  longitude 
associated  with  7.  Let  p , 0v,  and  \v  be  the  geocentric  distance,  geocentric  latitude, 
and  longitude  of  my.  Further,  let  be  the  angle  between  7 and  pv.  Now, 


cos  7 = sin  0 sin  6 +cosflcos0  cos  (X  - X ) (77) 

' V V V v V 


and  from  Page  115  of  Reference  14, 


Pn(cos7y)  = Pn(sin  0)Pn(sin  0^) 

+ 2 2 — — P"'  (sin  0)Pm  (sin  6 ) cos  m(X  - X ) 

m = l (n  + m)!  n n v 


K = ^ Jl()  ^rjP^sin^Pnfsin^) 
n (n  - m)! 

2 2 Pm  (sin  0)Pm  (sin  0 ) cos  m(X  - X 

m = iL(n  + m)!  n nr 


At  this  point,  we  separated  the  potential  tide  into  products  of  terms  which  depend 
on  r only  and  terms  which  contain  the  point  mass  coordinates  only: 


~ P„(sin  0)  « n 


terms  which  depend  terms  which  contain 

only  on  source  point  only  coordinates  of  the 
coordinates:  - general  point:  - 


From  here  results  the  following  sequence  of  definitions  and  equations. 


<>v  = tl  + <^2  +<^r3 

~ P (sin  0) 

<l>  . = H 2 pnP  (sin  6 ) ~ 

v 1 *v  n_  Q r\>  n V rn  + \ 


» " (n  - m)!  I P"1  (sin  0 ) I 

0 = 2p  2 2 -p"Pm  (sin  0 ) cos  mX  { — r; — cos  mX  / (84) 

v2  " n =o  m = i (n  + m)!  Hv  " v " r^ 


0 , = 2p  2 2 ~’pnPm(sin0  ) sin  mX  { n . — sin  mX  J (85) 

"3  v n =o  m = l (n  + m)!  v " v " rn  + 1 


P™  (sin  0) 


2 E”W 


* 


I 


0V2  = \ Fnn,Unn, 

vi  n = 0 in  = 1 


0 , = V 2 H^„,  V 


vi  n“o  m~l  "nm  ',m 


E"  = — — P (sin  0 ) 

n M Rn  n 


tlL  (l1  ~ m)!  £l  pn 


Fp  =2—  : — — F"(sin0  ) cos  mX 

/i  (n  + m)!  Rn  " 


= Mr  (n  Pm  (sin  9 ) sin  mX 

»"*  P (n  + m)!  Rn  n 


mR" 

wn  = prr,  pn<si"°> 


u 


AiRn 


n rn  rn  + ' n 


Pm  (sin  0)  cos  mX 


uR" 

V — — P™  (sin  0)  sin  mX 


n ni  j.n  ♦ 1 n 


= i;  } ' W + 2 v Fv  U + 22  H1’  V 

r fi  . 1 nm  nm  , nm  nm 


n = 0 


n =0  m = 1 


n = 0 m = 1 


Now.  because 


V „ = 0 
n 0 


US 


(sin  tnX)m  =n  = 0 


A -3 


i- ■ •'  ... 


•,!«>%  . A *. 


(87) 
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(93) 


(94) 


(95) 


(96) 

(97) 


A 


there  results 


v 2 H"  V 

n = 0 m = I nm  nm 


= 1' 


2 H"  V 

n = 0 m = 0 nm  nm 


(98) 


NovV,  redefine 


F*' 

n m 


F" 

nm 


= (2  - 6°  ) — (n~  m)!  — r (sin  0 ) cos  mX 
m ju  (n  + m)!  Rn  " 


(99) 


we  have 


2 E"  W +2  2 F"  U = 


n = 0 


n = 0 m = 1 


nm  nm 


^ £ F"  Unm 

n = 0 m = 0 nm  m 


(100) 


Finally,  remember  that 


0 = 2 0„ 

V-  1 


(101) 


and  also  consider  the  fact  that  whenever  we  extended  the  summation  over  n to 
infinity,  we  actually  thought  in  terms  of  a finite  sum,  extended  to  a finite  number, 
n . There  results  now  the  following  collection  of  equations,  for  the  tidal 

potential : 


max  n 

0=2  2 (F  U + H V ) 

v ' nm  nm  nm  nm' 

n = 0 m =0 


(102) 


= 2 

F” 

(103) 

v=  1 
N 

nm 

= 2 

H" 

(104) 

p=  l 

nm 
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(105) 


ir 

n m 


„ /a„  (n  - m)!  p'l 

- (2-  6 ) — — -Pm(sin0  ) cos  in X 

1,1  p (n  + m)!  R"  n 

P„  (n  - in)!  p" 

= — - P'"  (sin  0 ) sin  m A 

p (n  + in)!  R"  n ‘‘ 

pR" 

= -^rj-  P"  (sin  0)  cos  mX 
pR'1 

= nt ! P™  (sin  0)  sin  mX 


(106) 

(107) 

(108) 


Two  tilings  are  obvious  when  contemplating  Equations  105  and  106.  Firstly, 
according  to  the  Equations  15,  16,  17,  and  75,  it  appears  that  F"m  and  lF'm  can 
be  written  as  sums  of  two  terms  each  of  which  is  a cosine  or  sine  function  of  a 
time-dependent  argument,  multiplied  by  an  expression  which  depends  on  point  mass 
position.  The  same  is  true  for  F and  H . Secondly,  P and  H"  contain  the 
point  mass  coordinates  in  a form  which  suggests  that  from  each  a solid  harmonic  of 
the  latter  coordinates  should  be  factored  out.  This  will  enable  us  to  apply  the 
recurrence  relations  which  will  greatlv  facilitate  the  task  of  evaluating  the  F1’  and 
H''m . Both  aspects  of  the  matter  were  reflected  when  formulating  the  computer 
algorithm  for  the  tide  potential  as  listed  in  the  main  body  of  the  report. 

As  far  as  the  recurrence  relationships  are  concerned  which  we  just  mentioned 
and  which  are  also  invoked  on  several  occasions  in  the  main  body  of  the  report, 
these  reflect  the  familiar  property  of  the  solid  harmonics.  They  have  in  the  past 
found  frequent  use  in  our  computer  programs  for  satellite  geodesy,  in  particular  on 
those  occasions  when  we  faced  the  task  of  evaluating  the  complicated  expressions 
associated  with  the  earth’s  gravity  field.  As  they  are  not  novel,  they  exist  in  the 
form  of  working  notes  only  and  were  never  formally  documented. 
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APPENDIX  B 

NOTES  ON  THE  ALGORITHM  FOR  LUNAR  POSITION 


NOTES  ON  THE  ALGORITHM  EOR  LUNAR  POSITION 


Included  in  the  chapter  on  the  tidal  point  masses  is  an  algorithm  for  the  mean 
longitude  of  the  moon.  This  consists  of  Equations  18  through  25.  It  specifies  values 
of  lunar  mean  longitude,  at  00. 001'  UT,  for  each  day  during  the  calendar  years 
1975,  1976,  and  1977.  A similar  algorithm  appears  in  Reference  4 (see  chapter  THE 
MEAN  LONGITUDES  OE  SUN  AND  MOON  and  Appendix  B).  It  differs  from  our 
present  mean  longitude  routine  insofar  as  it  represents  the  mean  longitudes  for  the 
sun  as  well  as  for  the  moon,  at  arbitrary  instants  of  Universal  time.  Both  algorithms 
are  special  cases  of  a more  general  computer  routine  which  concerns  the  mean 
longitudes  of  sun,  moon,  and  lunar  perigee  and  which  was  formulated  by  one  of 
the  authors  (Groeger)  for  use  with  the  Schwiderski  Ocean  Tide  Program.  No  formal 
documentation  is  available  for  the  latter  routine.  But.  we  expect  to  publish  it 
sometime  in  the  future  in  a suitable  context. 

It  may  be  desirable  to  state  the  precise  meaning  of  the  term  “lunar  mean 
longitude.”  s (this  symbol  does  not  appear  elsewhere  in  this  report)  is  the  mean 
mean  longitude  of  the  moon.  Tin's  is  the  angle,  from  the  mean  equinox  of  date, 
measured  along  the  mean  ecliptic  of  date,  to  the  mean  ascending  node,  plus  the 
angle  measured  from  that  latter  point,  along  the  mean  orbit,  to  the  mean  perigee, 
plus  the  instantaneous  value  of  the  mean  mean  anomaly  corresponding  to  the 
position  of  the  moon,  s is  one  of  the  mean  elements  occurring  in  the  theory  of  the 
perturbed  motion,  relative  to  the  earth,  of  the  moon.  As  such,  it  is  the  secular  and 
the  very  long  periodic  portion  of  the  lunar  longitude.  The  latter  is  also  an  orbital 
element  in  lunar  theory,  x as  resulting  from  the  present  version  of  the  lunar  mean 
longitude  algorithm  is,  as  already  said  above,  equal  to  s.  the  latter  being  evaluated 
at  midnight  UT. 

Our  References  15.  16  (especially  Pages  537-540),  and  17  (Pages  98  and  107) 
refer  to  s as  a “mean”  longitude.  This  is  astronomical  usage,  adopted  for  the  sake 
of  brevity.  According  to  the  above  definition  of  s,  the  physical  item  associated  with 
s is  actually  a mean  mean  longitude. 

Further  elaborations  related  to  the  mean  longitude  algorithm  may  be  found  in 
Appendix  B of  Reference  4.  In  addition,  several  figures  will  be  found  there  which 
illustrate  the  concept  of  the  Julian  day,  for  the  benefit  of  those  who  only 
occasionally  encounter  that  unit  of  time. 
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Finally,  our  Equation  25  for  x has  an  accuracy  potential  of  several  hundred 
meters,  in  the  position  of  the  lunar  sub  point  on  the  earth’s  surface.  We  are  aware 
that  this  is  excessively  precise  for  the  purposes  of  present  tide  models.  It  was, 
however,  decided  to  render  the  polynomial  exactly  as  obtained  from  Reference  16, 
because  this  will  anticipate  possible  further  requirements.  We  do  not  expect  that  any 
worthwhile  economies  in  computer  time  would  result  if  a truncated  polynomial  was 
used. 
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